Import packages¶
In [1]:
import sys
print(sys.version)
import numpy as np
import matplotlib.pyplot as plt
import torch
from torch import nn
from rasterio.plot import show
from datetime import datetime, timedelta
import os
import pandas as pd
import rasterio
# Get today's date
today_date = datetime.today().strftime('%Y-%m-%d')
3.12.5 | packaged by Anaconda, Inc. | (main, Sep 12 2024, 18:18:29) [MSC v.1929 64 bit (AMD64)]
Import the tabular data¶
In [2]:
buffalo_id = 2005
n_samples = 10297 # 2005 has 10297 samples
# buffalo_id = 2018
# n_samples = 9440 # 2018 has 9440 samples
# Specify the path to your CSV file
csv_file_path = f'../buffalo_local_data_id/buffalo_{buffalo_id}_data_df_lag_1hr_n{n_samples}.csv'
# Read the CSV file into a DataFrame
buffalo_df = pd.read_csv(csv_file_path)
# Display the first few rows of the DataFrame
print(buffalo_df.head())
x_ y_ t_ id x1_ \
0 41969.310875 -1.435671e+06 2018-07-25T01:04:23Z 2005 41969.310875
1 41921.521939 -1.435654e+06 2018-07-25T02:04:39Z 2005 41921.521939
2 41779.439594 -1.435601e+06 2018-07-25T03:04:17Z 2005 41779.439594
3 41841.203272 -1.435635e+06 2018-07-25T04:04:39Z 2005 41841.203272
4 41655.463332 -1.435604e+06 2018-07-25T05:04:27Z 2005 41655.463332
y1_ x2_ y2_ x2_cent y2_cent ... \
0 -1.435671e+06 41921.521939 -1.435654e+06 -47.788936 16.857110 ...
1 -1.435654e+06 41779.439594 -1.435601e+06 -142.082345 53.568427 ...
2 -1.435601e+06 41841.203272 -1.435635e+06 61.763677 -34.322938 ...
3 -1.435635e+06 41655.463332 -1.435604e+06 -185.739939 31.003534 ...
4 -1.435604e+06 41618.651923 -1.435608e+06 -36.811409 -4.438037 ...
bearing bearing_sin bearing_cos ta cos_ta x_min \
0 2.802478 0.332652 -0.943050 1.367942 0.201466 40706.810875
1 2.781049 0.352783 -0.935705 -0.021429 0.999770 40659.021939
2 -0.507220 -0.485749 0.874098 2.994917 -0.989262 40516.939594
3 2.976198 0.164641 -0.986354 -2.799767 -0.942144 40578.703272
4 -3.021610 -0.119695 -0.992811 0.285377 0.959556 40392.963332
x_max y_min y_max ndvi_index
0 43231.810875 -1.436934e+06 -1.434409e+06 8
1 43184.021939 -1.436917e+06 -1.434392e+06 8
2 43041.939594 -1.436863e+06 -1.434338e+06 8
3 43103.703272 -1.436898e+06 -1.434373e+06 8
4 42917.963332 -1.436867e+06 -1.434342e+06 8
[5 rows x 32 columns]
Importing spatial data¶
Global layers¶
NDVI¶
In [3]:
# for monthly NDVI
file_path = '../mapping/cropped rasters/ndvi_monthly.tif'
# read the raster file
with rasterio.open(file_path) as src:
# Read the raster band as separate variable
ndvi_global = src.read([i for i in range(1, src.count + 1)])
# Get the metadata of the raster
ndvi_meta = src.meta
raster_transform = src.transform
# Print the metadata to check for time component
print("Metadata:", ndvi_meta)
# Check for specific time-related metadata
if 'TIFFTAG_DATETIME' in src.tags():
print("Time component found:", src.tags()['TIFFTAG_DATETIME'])
else:
print("No explicit time component found in metadata.")
# the rasters don't contain a time component
Metadata: {'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 24, 'crs': CRS.from_epsg(3112), 'transform': Affine(25.0, 0.0, 0.0,
0.0, -25.0, -1406000.0)}
No explicit time component found in metadata.
In [4]:
print(ndvi_meta)
print(raster_transform)
print(ndvi_global.shape)
# Replace NaNs in the original array with -1, which represents water
ndvi_global = np.nan_to_num(ndvi_global, nan=-1.0)
# from the stack of local layers
ndvi_max = 0.8220
ndvi_min = -0.2772
ndvi_global_tens = torch.from_numpy(ndvi_global)
# Normalizing the data
ndvi_global_norm = (ndvi_global_tens - ndvi_min) / (ndvi_max - ndvi_min)
# plt.imshow(ndvi_global_norm.numpy())
# plt.colorbar()
# plt.show()
plt.imshow(ndvi_global_norm[1,:,:].numpy())
plt.colorbar()
plt.show()
plt.imshow(ndvi_global_norm[8,:,:].numpy())
plt.colorbar()
plt.show()
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 24, 'crs': CRS.from_epsg(3112), 'transform': Affine(25.0, 0.0, 0.0,
0.0, -25.0, -1406000.0)}
| 25.00, 0.00, 0.00|
| 0.00,-25.00,-1406000.00|
| 0.00, 0.00, 1.00|
(24, 2280, 2400)
Canopy cover¶
In [5]:
file_path = '../mapping/cropped rasters/canopy_cover.tif'
# read the raster file
with rasterio.open(file_path) as src:
# Read the raster band as separate variable
canopy_global = src.read(1)
# Get the metadata of the raster
canopy_meta = src.meta
In [6]:
print(canopy_meta)
print(canopy_global.shape)
# from the stack of local layers
canopy_max = 82.5000
canopy_min = 0.0
canopy_global_tens = torch.from_numpy(canopy_global)
# Normalizing the data
canopy_global_norm = (canopy_global_tens - canopy_min) / (canopy_max - canopy_min)
plt.imshow(canopy_global_norm.numpy())
plt.colorbar()
plt.show()
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_epsg(3112), 'transform': Affine(25.0, 0.0, 0.0,
0.0, -25.0, -1406000.0)}
(2280, 2400)
Herbaceous vegetation¶
In [7]:
file_path = '../mapping/cropped rasters/veg_herby.tif'
# read the raster file
with rasterio.open(file_path) as src:
# Read the raster band as separate variable
herby_global = src.read(1)
# Get the metadata of the raster
herby_meta = src.meta
In [8]:
print(herby_meta)
print(herby_global.shape)
# from the stack of local layers
herby_max = 1.0
herby_min = 0.0
herby_global_tens = torch.from_numpy(herby_global)
# Normalizing the data
herby_global_norm = (herby_global_tens - herby_min) / (herby_max - herby_min)
plt.imshow(herby_global_norm.numpy())
plt.colorbar()
plt.show()
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.3999999521443642e+38, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_epsg(3112), 'transform': Affine(25.0, 0.0, 0.0,
0.0, -25.0, -1406000.0)}
(2280, 2400)
Slope¶
In [9]:
file_path = '../mapping/cropped rasters/slope.tif'
# read the raster file
with rasterio.open(file_path) as src:
# Read the raster band as separate variable
slope_global = src.read(1)
# Get the metadata of the raster
slope_meta = src.meta
slope_transform = src.transform # same as the raster transform in the NDVI raster read
print(slope_global)
print(slope_transform)
[[ nan nan nan ... nan nan nan] [ nan 0.3352837 0.39781624 ... 0. 0. nan] [ nan 0.3983888 0.48142004 ... 0. 0. nan] ... [ nan 2.215875 1.9798415 ... 1.5078747 1.268342 nan] [ nan 1.9740707 1.7354656 ... 1.697194 1.4880029 nan] [ nan nan nan ... nan nan nan]] | 25.00, 0.00, 0.00| | 0.00,-25.00,-1406000.00| | 0.00, 0.00, 1.00|
In [10]:
print(slope_meta)
print(slope_global.shape)
# Replace NaNs in the original array with -1, which represents water
slope_global = np.nan_to_num(slope_global, nan=0.0)
# from the stack of local layers
slope_max = 12.2981
slope_min = 0.0006
slope_global_tens = torch.from_numpy(slope_global)
# Normalizing the data
slope_global_norm = (slope_global_tens - slope_min) / (slope_max - slope_min)
plt.imshow(slope_global_norm.numpy())
plt.colorbar()
plt.show()
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': nan, 'width': 2400, 'height': 2280, 'count': 1, 'crs': CRS.from_epsg(3112), 'transform': Affine(25.0, 0.0, 0.0,
0.0, -25.0, -1406000.0)}
(2280, 2400)
Testing functions for selecting subsets of the raster layer, using torch objects
In [11]:
# Create a figure and axis with matplotlib
fig, ax = plt.subplots(figsize=(7.5, 7.5))
# Plot the raster
show(slope_global, transform=raster_transform, ax=ax, cmap='viridis')
# Set the title and labels
ax.set_title('Raster with Geographic Coordinates')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Show the plot
plt.show()
In [12]:
x, y = 5.9e4, -1.447e6
print(x, y)
# Convert geographic coordinates to pixel coordinates
px, py = ~raster_transform * (x, y)
# Round pixel coordinates to integers
px, py = int(round(px)), int(round(py))
# Print the pixel coordinates
print(px, py)
59000.0 -1447000.0 2360 1640
In [13]:
window_size = 101
# Define half the window size
half_window = window_size // 2
# Calculate the window boundaries
row_start = py - half_window
row_stop = py + half_window + 1
col_start = px - half_window
col_stop = px + half_window + 1
# Initialize the subset array with zeros (or any other padding value)
subset = np.zeros((window_size, window_size), dtype=slope_global.dtype)
# Calculate the valid region within the raster bounds
valid_row_start = max(0, row_start)
valid_row_stop = min(slope_global.shape[0], row_stop)
valid_col_start = max(0, col_start)
valid_col_stop = min(slope_global.shape[1], col_stop)
# Calculate the corresponding region in the subset array
subset_row_start = valid_row_start - row_start
subset_row_stop = subset_row_start + (valid_row_stop - valid_row_start)
subset_col_start = valid_col_start - col_start
subset_col_stop = subset_col_start + (valid_col_stop - valid_col_start)
# Copy the valid region from the raster array to the subset array
subset[subset_row_start:subset_row_stop, subset_col_start:subset_col_stop] = \
slope_global[valid_row_start:valid_row_stop, valid_col_start:valid_col_stop]
Plot the raster layer
In [14]:
# plot the subset
plt.imshow(subset, cmap='viridis')
plt.colorbar()
plt.show()
Subset function¶
In [15]:
def subset_raster_with_padding_torch(raster_tensor, x, y, window_size, transform):
# Convert geographic coordinates to pixel coordinates
px, py = ~transform * (x, y)
# Round pixel coordinates to integers
px, py = int(round(px)), int(round(py))
# Define half the window size
half_window = window_size // 2
# Calculate the window boundaries
row_start = py - half_window
row_stop = py + half_window + 1
col_start = px - half_window
col_stop = px + half_window + 1
# Initialize the subset tensor with zeros (or any other padding value)
subset = torch.full((window_size, window_size), -1.0, dtype=raster_tensor.dtype)
# Calculate the valid region within the raster bounds
valid_row_start = max(0, row_start)
valid_row_stop = min(raster_tensor.shape[0], row_stop)
valid_col_start = max(0, col_start)
valid_col_stop = min(raster_tensor.shape[1], col_stop)
# Calculate the corresponding region in the subset tensor
subset_row_start = valid_row_start - row_start
subset_row_stop = subset_row_start + (valid_row_stop - valid_row_start)
subset_col_start = valid_col_start - col_start
subset_col_stop = subset_col_start + (valid_col_stop - valid_col_start)
# Copy the valid region from the raster tensor to the subset tensor
subset[subset_row_start:subset_row_stop, subset_col_start:subset_col_stop] = \
raster_tensor[valid_row_start:valid_row_stop, valid_col_start:valid_col_stop]
return subset, col_start, row_start
Testing the subset function
In [16]:
x = 5.9e4
y = -1.41e6
window_size = 101
which_ndvi = 1
ndvi_subset, origin_x, origin_y = subset_raster_with_padding_torch(ndvi_global_norm[which_ndvi,:,:], x, y, window_size, raster_transform)
canopy_subset, origin_x, origin_y = subset_raster_with_padding_torch(canopy_global_norm, x, y, window_size, raster_transform)
herby_subset, origin_x, origin_y = subset_raster_with_padding_torch(herby_global_norm, x, y, window_size, raster_transform)
slope_subset, origin_x, origin_y = subset_raster_with_padding_torch(slope_global_norm, x, y, window_size, raster_transform)
# Plot the subset
fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(ndvi_subset.numpy(), cmap='viridis')
axs[0, 0].set_title('NDVI Subset')
# axs[0, 0].axis('off')
axs[0, 1].imshow(canopy_subset.numpy(), cmap='viridis')
axs[0, 1].set_title('Canopy Cover Subset')
# axs[0, 1].axis('off')
axs[1, 0].imshow(herby_subset.numpy(), cmap='viridis')
axs[1, 0].set_title('Herbaceous Vegetation Subset')
# axs[1, 0].axis('off')
axs[1, 1].imshow(slope_subset.numpy(), cmap='viridis')
axs[1, 1].set_title('Slope Subset')
# axs[1, 1].axis('off')
Out[16]:
Text(0.5, 1.0, 'Slope Subset')
In [17]:
x = 3.7e4
y = -1.435e6
print(x, y)
# Convert geographic coordinates to pixel coordinates
px, py = ~raster_transform * (x, y)
# Round pixel coordinates to integers
px, py = int(round(px)), int(round(py))
# Print the pixel coordinates
print(px, py)
37000.0 -1435000.0 1480 1160
Running the model on the subset layers¶
Set the device for the model¶
In [18]:
# train on the GPU or on the CPU, if a GPU is not available
device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')
print(f"Using {device} device")
Using cpu device
Define the model¶
In [19]:
class Conv2d_block_toFC(nn.Module):
def __init__(self, params):
super(Conv2d_block_toFC, self).__init__()
self.batch_size = params.batch_size
self.input_channels = params.input_channels
self.output_channels = params.output_channels
self.kernel_size = params.kernel_size
self.stride = params.stride
self.kernel_size_mp = params.kernel_size_mp
self.stride_mp = params.stride_mp
self.padding = params.padding
self.image_dim = params.image_dim
self.device = params.device
self.conv2d = nn.Sequential(
nn.Conv2d(in_channels=self.input_channels, out_channels=self.output_channels, kernel_size=self.kernel_size, stride=self.stride, padding=self.padding),
nn.ReLU(),
nn.MaxPool2d(kernel_size=self.kernel_size_mp, stride=self.stride_mp),
nn.Conv2d(in_channels=self.output_channels, out_channels=self.output_channels, kernel_size=self.kernel_size, stride=self.stride, padding=self.padding),
nn.ReLU(),
nn.MaxPool2d(kernel_size=self.kernel_size_mp, stride=self.stride_mp),
nn.Flatten())
def forward(self, x):
return self.conv2d(x)
class Conv2d_block_spatial(nn.Module):
def __init__(self, params):
super(Conv2d_block_spatial, self).__init__()
self.batch_size = params.batch_size
self.input_channels = params.input_channels
self.output_channels = params.output_channels
self.kernel_size = params.kernel_size
self.stride = params.stride
# self.kernel_size_mp = params.kernel_size_mp
# self.stride_mp = params.stride_mp
self.padding = params.padding
self.image_dim = params.image_dim
self.device = params.device
self.conv2d = nn.Sequential(
nn.Conv2d(in_channels=self.input_channels, out_channels=self.output_channels, kernel_size=self.kernel_size, stride=self.stride, padding=self.padding),
nn.ReLU(),
nn.Conv2d(in_channels=self.output_channels, out_channels=self.output_channels, kernel_size=self.kernel_size, stride=self.stride, padding=self.padding),
nn.ReLU(),
nn.Conv2d(in_channels=self.output_channels, out_channels=1, kernel_size=self.kernel_size, stride=self.stride, padding=self.padding)
)
def forward(self, x):
conv2d_spatial = self.conv2d(x).squeeze(dim = 1)
conv2d_spatial = conv2d_spatial - torch.logsumexp(conv2d_spatial, dim = (1, 2), keepdim = True)
return conv2d_spatial
class FCN_block_all_habitat(nn.Module):
def __init__(self, params):
super(FCN_block_all_habitat, self).__init__()
self.batch_size = params.batch_size
self.dense_dim_in_all = params.dense_dim_in_all
self.dense_dim_hidden = params.dense_dim_hidden
self.dense_dim_out = params.dense_dim_out
self.image_dim = params.image_dim
self.device = params.device
self.dropout = params.dropout
self.ffn = nn.Sequential(
nn.Linear(self.dense_dim_in_all, self.dense_dim_hidden),
nn.Dropout(self.dropout),
nn.ReLU(),
nn.Linear(self.dense_dim_hidden, self.dense_dim_hidden),
nn.Dropout(self.dropout),
nn.ReLU(),
nn.Linear(self.dense_dim_hidden, self.image_dim * self.image_dim)
)
def forward(self, x):
return self.ffn(x)
class FCN_block_all_movement(nn.Module):
def __init__(self, params):
super(FCN_block_all_movement, self).__init__()
self.batch_size = params.batch_size
self.dense_dim_in_all = params.dense_dim_in_all
self.dense_dim_hidden = params.dense_dim_hidden
self.dense_dim_out = params.dense_dim_out
self.image_dim = params.image_dim
self.device = params.device
self.num_movement_params = params.num_movement_params
self.dropout = params.dropout
self.ffn = nn.Sequential(
nn.Linear(self.dense_dim_in_all, self.dense_dim_hidden),
nn.Dropout(self.dropout),
nn.ReLU(),
nn.Linear(self.dense_dim_hidden, self.dense_dim_hidden),
nn.Dropout(self.dropout),
nn.ReLU(),
nn.Linear(self.dense_dim_hidden, self.num_movement_params)
)
def forward(self, x):
return self.ffn(x)
##################################################
# Mixture of two Gamma and von Mises distributions with the von Mises mu parameters allowed to vary
##################################################
class Params_to_Grid_Block(nn.Module):
def __init__(self, params):
super(Params_to_Grid_Block, self).__init__()
self.batch_size = params.batch_size
self.image_dim = params.image_dim
self.pixel_size = params.pixel_size
self.center = self.image_dim // 2
y, x = np.indices((self.image_dim, self.image_dim))
self.distance_layer = torch.from_numpy(np.sqrt((self.pixel_size*(x - self.center))**2 + (self.pixel_size*(y - self.center))**2))
# change the centre cell to the average distance from the centre to the edge of the pixel
self.distance_layer[self.center, self.center] = 0.56*self.pixel_size # average distance from the centre to the perimeter of the pixel (accounting for longer distances at the corners)
# self.bearing_layer = torch.from_numpy(np.arctan2(y - self.center, x - self.center))
self.bearing_layer = torch.from_numpy(np.arctan2(self.center - y, x - self.center))
self.device = params.device
# Gamma desnities for the mixture distribution
def gamma_density(self, x, shape, scale):
# Ensure all tensors are on the same device as x
shape = shape.to(x.device)
scale = scale.to(x.device)
return -1*torch.lgamma(shape) -shape*torch.log(scale) + (shape - 1)*torch.log(x) - x/scale
# von Mises densities for the mixture distribution
def vonmises_density(self, x, kappa, vm_mu):
# Ensure all tensors are on the same device as x
kappa = kappa.to(x.device)
vm_mu = vm_mu.to(x.device)
return kappa*torch.cos(x - vm_mu) - 1*(np.log(2*torch.pi) + torch.log(torch.special.i0(kappa)))
def forward(self, x, bearing):
# parameters of the first mixture distribution
gamma_shape1 = torch.exp(x[:, 0]).unsqueeze(0).unsqueeze(0)
gamma_shape1 = gamma_shape1.repeat(self.image_dim, self.image_dim, 1)
gamma_shape1 = gamma_shape1.permute(2, 0, 1)
gamma_scale1 = torch.exp(x[:, 1]).unsqueeze(0).unsqueeze(0)
gamma_scale1 = gamma_scale1.repeat(self.image_dim, self.image_dim, 1)
gamma_scale1 = gamma_scale1.permute(2, 0, 1)
gamma_weight1 = torch.exp(x[:, 2]).unsqueeze(0).unsqueeze(0)
gamma_weight1 = gamma_weight1.repeat(self.image_dim, self.image_dim, 1)
gamma_weight1 = gamma_weight1.permute(2, 0, 1)
# parameters of the second mixture distribution
gamma_shape2 = torch.exp(x[:, 3]).unsqueeze(0).unsqueeze(0)
gamma_shape2 = gamma_shape2.repeat(self.image_dim, self.image_dim, 1)
gamma_shape2 = gamma_shape2.permute(2, 0, 1)
gamma_scale2 = torch.exp(x[:, 4]).unsqueeze(0).unsqueeze(0)
gamma_scale2 = gamma_scale2 * 500 ### to transform the scale parameter to be near 1
gamma_scale2 = gamma_scale2.repeat(self.image_dim, self.image_dim, 1)
gamma_scale2 = gamma_scale2.permute(2, 0, 1)
gamma_weight2 = torch.exp(x[:, 5]).unsqueeze(0).unsqueeze(0)
gamma_weight2 = gamma_weight2.repeat(self.image_dim, self.image_dim, 1)
gamma_weight2 = gamma_weight2.permute(2, 0, 1)
# Apply softmax to the weights
gamma_weights = torch.stack([gamma_weight1, gamma_weight2], dim=0)
gamma_weights = torch.nn.functional.softmax(gamma_weights, dim=0)
gamma_weight1 = gamma_weights[0]
gamma_weight2 = gamma_weights[1]
# calculation of Gamma densities
gamma_density_layer1 = self.gamma_density(self.distance_layer, gamma_shape1, gamma_scale1).to(device)
gamma_density_layer2 = self.gamma_density(self.distance_layer, gamma_shape2, gamma_scale2).to(device)
# combining both densities to create a mixture distribution using logsumexp
logsumexp_gamma_corr = torch.max(gamma_density_layer1, gamma_density_layer2)
gamma_density_layer = logsumexp_gamma_corr + torch.log(gamma_weight1 * torch.exp(gamma_density_layer1 - logsumexp_gamma_corr) + gamma_weight2 * torch.exp(gamma_density_layer2 - logsumexp_gamma_corr))
## Von Mises Distributions
# calculate the new bearing from the turning angle
# takes in the bearing from the previous step and adds the turning angle
bearing_new1 = x[:, 6] + bearing[:, 0]
# the new bearing becomes the mean of the von Mises distribution
# the estimated parameter [x:, 6] is the turning angle of the next step
# which is always in reference to the input bearing
vonmises_mu1 = bearing_new1.unsqueeze(0).unsqueeze(0)
vonmises_mu1 = vonmises_mu1.repeat(self.image_dim, self.image_dim, 1)
vonmises_mu1 = vonmises_mu1.permute(2, 0, 1)
# parameters of the first von Mises distribution
vonmises_kappa1 = torch.exp(x[:, 7]).unsqueeze(0).unsqueeze(0)
vonmises_kappa1 = vonmises_kappa1.repeat(self.image_dim, self.image_dim, 1)
vonmises_kappa1 = vonmises_kappa1.permute(2, 0, 1)
vonmises_weight1 = torch.exp(x[:, 8]).unsqueeze(0).unsqueeze(0)
vonmises_weight1 = vonmises_weight1.repeat(self.image_dim, self.image_dim, 1)
vonmises_weight1 = vonmises_weight1.permute(2, 0, 1)
# vm_mu and weight for the second von Mises distribution
bearing_new2 = x[:, 9] + bearing[:, 0]
vonmises_mu2 = bearing_new2.unsqueeze(0).unsqueeze(0)
vonmises_mu2 = vonmises_mu2.repeat(self.image_dim, self.image_dim, 1)
vonmises_mu2 = vonmises_mu2.permute(2, 0, 1)
# parameters of the second von Mises distribution
vonmises_kappa2 = torch.exp(x[:, 10]).unsqueeze(0).unsqueeze(0)
vonmises_kappa2 = vonmises_kappa2.repeat(self.image_dim, self.image_dim, 1)
vonmises_kappa2 = vonmises_kappa2.permute(2, 0, 1)
vonmises_weight2 = torch.exp(x[:, 11]).unsqueeze(0).unsqueeze(0)
vonmises_weight2 = vonmises_weight2.repeat(self.image_dim, self.image_dim, 1)
vonmises_weight2 = vonmises_weight2.permute(2, 0, 1)
# Apply softmax to the weights
vonmises_weights = torch.stack([vonmises_weight1, vonmises_weight2], dim=0)
vonmises_weights = torch.nn.functional.softmax(vonmises_weights, dim=0)
vonmises_weight1 = vonmises_weights[0]
vonmises_weight2 = vonmises_weights[1]
# calculation of von Mises densities
vonmises_density_layer1 = self.vonmises_density(self.bearing_layer, vonmises_kappa1, vonmises_mu1).to(device)
vonmises_density_layer2 = self.vonmises_density(self.bearing_layer, vonmises_kappa2, vonmises_mu2).to(device)
# combining both densities to create a mixture distribution using the logsumexp trick
logsumexp_vm_corr = torch.max(vonmises_density_layer1, vonmises_density_layer2)
vonmises_density_layer = logsumexp_vm_corr + torch.log(vonmises_weight1 * torch.exp(vonmises_density_layer1 - logsumexp_vm_corr) + vonmises_weight2 * torch.exp(vonmises_density_layer2 - logsumexp_vm_corr))
# combining the two distributions
movement_grid = gamma_density_layer + vonmises_density_layer # Gamma and von Mises densities are on the log-scale
# normalise the movement predictions
movement_grid = movement_grid - torch.logsumexp(movement_grid, dim = (1, 2), keepdim = True)
return movement_grid
class Scalar_to_Grid_Block(nn.Module):
def __init__(self, params):
super(Scalar_to_Grid_Block, self).__init__()
self.batch_size = params.batch_size
self.image_dim = params.image_dim
self.device = params.device
def forward(self, x):
num_scalars = x.shape[1]
scalar_map = x.view(x.shape[0], num_scalars, 1, 1).expand(x.shape[0], num_scalars, self.image_dim, self.image_dim)
return scalar_map
class ConvJointModel(nn.Module):
def __init__(self, params):
super(ConvJointModel, self).__init__()
self.scalar_grid_output = Scalar_to_Grid_Block(params)
self.conv_habitat = Conv2d_block_spatial(params)
self.conv_movement = Conv2d_block_toFC(params)
self.fcn_movement_all = FCN_block_all_movement(params)
self.movement_grid_output = Params_to_Grid_Block(params)
self.device = params.device
def forward(self, x):
spatial_data_x = x[0]
scalars_to_grid = x[1]
bearing_x = x[2]
# SCALAR GRIDS
scalar_grids = self.scalar_grid_output(scalars_to_grid)
all_spatial = torch.cat([spatial_data_x, scalar_grids], dim = 1)
# print(f"Shape after scalar grid: {all_spatial.shape}") # Debugging print
# HABITAT SELECTION
output_habitat = self.conv_habitat(all_spatial)
# print(f"Shape after CNN habitat: {output_habitat.shape}") # Debugging print
# MOVEMENT
conv_movement = self.conv_movement(all_spatial)
# print(f"Shape after CNN to FC movement: {conv_movement.shape}") # Debugging print
output_movement = self.fcn_movement_all(conv_movement)
# print(f"Shape after fcn_movement_all: {output_movement.shape}") # Debugging print
output_movement = self.movement_grid_output(output_movement, bearing_x)
# print(f"Shape after CNN movement: {output_movement.shape}") # Debugging print
# combine the habitat and movement predictions
output = torch.stack((output_habitat, output_movement), dim = -1)
return output
class ModelParams():
def __init__(self, dict_params):
self.batch_size = dict_params["batch_size"]
self.image_dim = dict_params["image_dim"]
self.pixel_size = dict_params["pixel_size"]
self.batch_size = dict_params["batch_size"]
self.dim_in_nonspatial_to_grid = dict_params["dim_in_nonspatial_to_grid"]
self.dense_dim_in_nonspatial = dict_params["dense_dim_in_nonspatial"]
self.dense_dim_hidden = dict_params["dense_dim_hidden"]
self.dense_dim_out = dict_params["dense_dim_out"]
self.batch_size = dict_params["batch_size"]
self.dense_dim_in_all = dict_params["dense_dim_in_all"]
self.dense_dim_hidden = dict_params["dense_dim_hidden"]
self.dense_dim_out = dict_params["dense_dim_out"]
self.batch_size = dict_params["batch_size"]
self.input_channels = dict_params["input_channels"]
self.output_channels = dict_params["output_channels"]
self.kernel_size = dict_params["kernel_size"]
self.stride = dict_params["stride"]
self.kernel_size_mp = dict_params["kernel_size_mp"]
self.stride_mp = dict_params["stride_mp"]
self.padding = dict_params["padding"]
self.image_dim = dict_params["image_dim"]
self.num_movement_params = dict_params["num_movement_params"]
self.dropout = dict_params["dropout"]
self.device = dict_params["device"]
Instantiate the model¶
In [20]:
params_dict = {"batch_size": 32,
"image_dim": 101, #number of pixels along the edge of each local patch/image
"pixel_size": 25, #number of metres along the edge of a pixel
"dim_in_nonspatial_to_grid": 4, #the number of scalar predictors that are converted to a grid and appended to the spatial features
"dense_dim_in_nonspatial": 4, #change this to however many other scalar predictors you have (bearing, velocity etc)
"dense_dim_hidden": 128, #number of nodes in the hidden layers
"dense_dim_out": 128, #number of nodes in the output of the fully connected block (FCN)
"dense_dim_in_all": 2500,# + 128, #number of inputs entering the fully connected block once the nonspatial features have been concatenated to the spatial features
"input_channels": 4 + 4, #number of spatial layers in each image + number of scalar layers that are converted to a grid
"output_channels": 4, #number of filters to learn
"kernel_size": 3, #the size of the 2D moving windows / kernels that are being learned
"stride": 1, #the stride used when applying the kernel. This reduces the dimension of the output if set to greater than 1
"kernel_size_mp": 2, #the size of the kernel that is used in max pooling operations
"stride_mp": 2, #the stride that is used in max pooling operations
"padding": 1, #the amount of padding to apply to images prior to applying the 2D convolution
"num_movement_params": 12, #number of parameters used to parameterise the movement kernel
"dropout": 0.1,
"device": device
}
params = ModelParams(params_dict)
model = ConvJointModel(params).to(device)
# print(model)
Load the model structure¶
In [49]:
# date of the trained model checkpoint
date = '2024-09-29'
# # load the model weights
# print(model.state_dict())
model.load_state_dict(torch.load(f'model_checkpoints/checkpoint_CNN_global_buffalo{buffalo_id}_TAmix_bearing-rev_adj_hab-covs_grid-only_{date}.pt', map_location=torch.device('cpu')))
# print(model.state_dict())
# model.eval()
C:\Users\for329\AppData\Local\Temp\ipykernel_24996\3300322696.py:6: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.
model.load_state_dict(torch.load(f'model_checkpoints/checkpoint_CNN_global_buffalo{buffalo_id}_TAmix_bearing-rev_adj_hab-covs_grid-only_{date}.pt', map_location=torch.device('cpu')))
Out[49]:
<All keys matched successfully>
Setup simulation parameters¶
Testing the subset function
In [50]:
# Create a mapping from day of the year to month index
def day_to_month_index(day_of_year):
# Calculate the year and the day within that year
base_date = datetime(2018, 1, 1)
date = base_date + timedelta(days=int(day_of_year) - 1)
year_diff = date.year - base_date.year
month_index = (date.month - 1) + (year_diff * 12) # month index (0-based, accounting for year change)
return month_index
yday = 60
month_index = day_to_month_index(yday)
print(month_index)
2
In [51]:
window_size = 101
# starting location of buffalo 2005
x = 41969.310875
y = -1.435671e+06
yday = 280
month_index = day_to_month_index(yday)
ndvi_subset, origin_x, origin_y = subset_raster_with_padding_torch(ndvi_global_norm[month_index,:,:], x, y, window_size, raster_transform)
canopy_subset, origin_x, origin_y = subset_raster_with_padding_torch(canopy_global_norm, x, y, window_size, raster_transform)
herby_subset, origin_x, origin_y = subset_raster_with_padding_torch(herby_global_norm, x, y, window_size, raster_transform)
slope_subset, origin_x, origin_y = subset_raster_with_padding_torch(slope_global_norm, x, y, window_size, raster_transform)
# Plot the subset
fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(ndvi_subset.numpy(), cmap='viridis')
axs[0, 0].set_title('NDVI')
# axs[0, 0].axis('off')
axs[0, 1].imshow(canopy_subset.numpy(), cmap='viridis')
axs[0, 1].set_title('Canopy Cover')
# axs[0, 1].axis('off')
axs[1, 0].imshow(herby_subset.numpy(), cmap='viridis')
axs[1, 0].set_title('Herbaceous Vegetation')
# axs[1, 0].axis('off')
axs[1, 1].imshow(slope_subset.numpy(), cmap='viridis')
axs[1, 1].set_title('Slope')
# axs[1, 1].axis('off')
Out[51]:
Text(0.5, 1.0, 'Slope')
Dealing with NaN values¶
In [52]:
layer = ndvi_subset
# layer = ndvi_global_norm
layer = slope_subset
# are there nans in the array?
# print(torch.isnan(layer).any())
# Create mask where original array has values of -1, which is only at the edges as everything else is normalised between 0 and 1
mask = np.where(layer == -1, -np.inf, 1)
plt.imshow(mask)
plt.colorbar()
plt.show()
# # Replace NaNs in the original array with 0
# layer_filled = np.nan_to_num(layer, nan=0.)
layer_updated = layer * mask
plt.imshow(layer_updated)
plt.colorbar()
plt.show()
Stack the subset layers¶
In [53]:
# Stack the channels along a new axis; here, 1 is commonly used for channel axis in PyTorch
subset_stack = torch.stack([ndvi_subset, canopy_subset, herby_subset, slope_subset], dim=0)
subset_stack
x1 = subset_stack
x1 = x1.unsqueeze(0)
print(x1.shape)
# print(x1)
torch.Size([1, 4, 101, 101])
Additional data¶
In [54]:
n_steps = 1000
def repeat_sequence(sequence, length_out):
return np.resize(sequence, length_out)
# hour of the day (hour) sequence
hour_t2 = np.resize(range(24), n_steps)
# print(hour_t2)
# convert to sine and cosine
hour_t2_sin = np.sin(2*np.pi*hour_t2/24)
hour_t2_cos = np.cos(2*np.pi*hour_t2/24)
# day of the year (yday) sequence
yday_sequence = np.repeat(range(yday, 365), 24)
yday_t2 = np.resize(yday_sequence, n_steps)
# print(yday_t2)
# convert to sine and cosine
yday_t2_sin = np.sin(2*np.pi*yday_t2/365)
yday_t2_cos = np.cos(2*np.pi*yday_t2/365)
# bearing vector
bearing = np.repeat(0, n_steps)
# Convert lists to PyTorch tensors
hour_t2_tensor = torch.tensor(hour_t2).float()
hour_t2_sin_tensor = torch.tensor(hour_t2_sin).float()
hour_t2_cos_tensor = torch.tensor(hour_t2_cos).float()
yday_t2_tensor = torch.tensor(yday_t2).float()
yday_t2_sin_tensor = torch.tensor(yday_t2_sin).float()
yday_t2_cos_tensor = torch.tensor(yday_t2_cos).float()
bearing = torch.tensor(bearing).float()
# Stack tensors column-wise
x2_full = torch.stack((hour_t2_sin_tensor, hour_t2_cos_tensor, yday_t2_sin_tensor, yday_t2_cos_tensor), dim=1)
# print(x2)
print(x2_full.shape)
# print(x2)
print(x2_full[59,:])
torch.Size([1000, 4]) tensor([ 0.2588, -0.9659, -0.9899, 0.1415])
Functions to recover the hour and yday from the sin and cosine terms¶
In [55]:
def recover_hour(sin_term, cos_term):
# Calculate the angle theta
theta = np.arctan2(sin_term, cos_term)
# Calculate hour_t2
hour = (12 * theta) / np.pi % 24
return hour
def recover_yday(sin_term, cos_term):
# Calculate the angle theta
theta = np.arctan2(sin_term, cos_term)
# Calculate hour_t2
yday = (365 * theta) / (2 * np.pi) % 365
return yday
Test the model on the subsetted layers¶
In [56]:
step = 20
test = model((x1, x2_full[step,:].unsqueeze(0), bearing[step].unsqueeze(0).unsqueeze(0)))
print(x2_full[step,:])
print(bearing[step].unsqueeze(0).unsqueeze(0))
x2_step = x2_full[step,:]
# Pull out the scalars
hour_t2_sin = x2_step.detach().numpy()[0]
hour_t2_cos = x2_step.detach().numpy()[1]
yday_t2_sin = x2_step.detach().numpy()[2]
yday_t2_cos = x2_step.detach().numpy()[3]
# recover the hour value for the step
hour_t2 = recover_hour(hour_t2_sin, hour_t2_cos)
print(hour_t2)
# recover the yday value for the step
yday_t2 = recover_yday(yday_t2_sin, yday_t2_cos)
print(yday_t2)
# plot the results of the habitat density as an image
hab_density = test.detach().numpy()[0,:,:,0]
plt.imshow(hab_density)
plt.colorbar()
plt.show()
# Create the mask for x and y coordinates
x_mask = np.ones_like(hab_density)
y_mask = np.ones_like(hab_density)
# mask out cells on the edges that affect the colour scale
x_mask[:, :3] = -np.inf
x_mask[:, 98:] = -np.inf
y_mask[:3, :] = -np.inf
y_mask[98:, :] = -np.inf
hab_density_mask = hab_density * x_mask * y_mask
plt.imshow(hab_density_mask)
plt.colorbar()
plt.show()
plt.imshow(np.exp(hab_density))
plt.colorbar()
plt.show()
# plot the results of the movement density as an image
plt.imshow(test.detach().numpy()[0,:,:,1])
plt.colorbar()
plt.show()
# plot the results of the exp movement density as an image
plt.imshow(np.exp(test.detach().numpy()[0,:,:,1]))
plt.colorbar()
plt.show()
test_cat = (test[:, :, :, 0] + test[:, :, :, 1])
test_cat = test_cat.squeeze()
print(test_cat.shape)
# # Create mask where original array has values of -1, which is only at the edges as everything else is normalised between 0 and 1
layer = ndvi_subset
mask = torch.where(layer == -1, torch.tensor(float('nan')), 1)
plt.imshow(mask)
plt.colorbar()
plt.show()
test_cat = test_cat * mask
test_cat = test_cat.detach().numpy()[:,:]
test_cat = test_cat * x_mask * y_mask
plt.imshow(test_cat)
plt.colorbar()
plt.show()
test_cat_exp = np.exp(test_cat)
test_cat_exp = test_cat_exp * x_mask * y_mask
plt.imshow(test_cat_exp)
plt.colorbar()
plt.show()
tensor([-0.8660, 0.5000, -0.9942, 0.1074]) tensor([[0.]]) 19.99999988868986 280.00000313675633
torch.Size([101, 101])
C:\Users\for329\AppData\Local\Temp\ipykernel_24996\4292268014.py:81: RuntimeWarning: invalid value encountered in multiply test_cat_exp = test_cat_exp * x_mask * y_mask
Sample from the probability surface¶
In [57]:
test_cat = (test[:, :, :, 0] + test[:, :, :, 1])
test_cat = test_cat.squeeze()
# sample from the array values
test_cat_exp = torch.exp(test_cat)
# replace NaNs with 0
test_cat_exp = torch.nan_to_num(test_cat_exp, nan=0.)
# normalise the array
test_cat_exp = test_cat_exp/torch.sum(test_cat_exp)
# print(test_cat_exp)
print(torch.sum(test_cat_exp))
# Flatten the probability surface
flat_prob_surface = test_cat_exp.flatten().detach().numpy()
print(flat_prob_surface)
# Generate the corresponding indices for the flattened array
indices = np.arange(flat_prob_surface.size)
print(indices)
# Sample from the flattened probability surface
sampled_index = np.random.choice(indices, p=flat_prob_surface)
print(sampled_index)
print(flat_prob_surface[sampled_index])
# Convert the sampled index back to 2D coordinates
sampled_coordinates = np.unravel_index(sampled_index, test_cat_exp.shape)
print("Sampled coordinates:", sampled_coordinates)
step_prob = np.exp(test_cat.detach().numpy()[:,:])
# Set the pixel of the next step, which is at (x, y) to 0
step_prob[sampled_coordinates] = -0.0000001
# step_prob[10, 30] = 1.0 # the y axis comes FIRST
plt.imshow(step_prob)
plt.colorbar()
plt.show()
tensor(1., dtype=torch.float64, grad_fn=<SumBackward0>) [5.29215473e-08 8.41615129e-08 1.32588328e-07 ... 4.26569001e-08 3.66841909e-08 4.81882274e-08] [ 0 1 2 ... 10198 10199 10200] 5304 0.0022398973703402613 Sampled coordinates: (52, 52)
Return sampled point to geographic coordinates¶
In [58]:
# original locations
print(x ,y)
print(sampled_coordinates)
print(px, py)
# row_start = py - half_window
# row_stop = py + half_window + 1
# col_start = px - half_window
# col_stop = px + half_window + 1
print(row_start, col_start)
# new_px = origin_x + sampled_coordinates[0]
# new_py = origin_y + sampled_coordinates[1]
# THE Y COORDINATE COMES FIRST
new_px = origin_x + sampled_coordinates[1]
new_py = origin_y + sampled_coordinates[0]
print(new_px)
print(new_py)
# Convert geographic coordinates to pixel coordinates
new_x, new_y = raster_transform * (new_px, new_py)
print(new_x, new_y)
41969.310875 -1435671.0 (52, 52) 45 58 1590 2310 1681 1189 42025.0 -1435725.0
Full trajectory function¶
In [59]:
def simulate_trajectory(global_raster_tensors,
scalars_to_grid,
# additional_data_tensor,
bearing,
window_size,
x_loc,
y_loc,
global_raster_transform):
results = [subset_raster_with_padding_torch(raster_tensor, x=x_loc, y=y_loc, window_size=window_size, transform=global_raster_transform) for raster_tensor in global_raster_tensors]
# Unpacking the results
subset_rasters_tensors, origin_xs, origin_ys = zip(*results)
# print(subset_rasters_tensors)
# Stack the processed tensors along a new dimension (e.g., dimension 0)
x1 = torch.stack(subset_rasters_tensors, dim=0)
x1 = x1.unsqueeze(0)
# print(x1.shape)
single_layer = x1[0, 0, :, :]
# plt.imshow(single_layer.detach().numpy()[:,:])
# plt.colorbar()
# plt.show()
# create masking layer to remove outside of the extent
mask = torch.where(single_layer == -1, torch.tensor(float('nan')), 1)
# test_cat_exp = torch.exp(test_cat)
# extract NaNs to pad and make masking layer
x2 = scalars_to_grid
# x3 = additional_data_tensor
x3 = bearing
test = model((x1, x2, x3))
# print(test.shape)
hab_log_prob = test[:, :, :, 0]
move_log_prob = test[:, :, :, 1]
step_log_prob = (hab_log_prob + move_log_prob)
step_log_prob = step_log_prob * mask
hab_log_prob = hab_log_prob.squeeze()
move_log_prob = move_log_prob.squeeze()
step_log_prob = step_log_prob.squeeze()
# sample from the array values
step_prob = torch.exp(step_log_prob)
step_prob = torch.nan_to_num(step_prob, nan=0.)
step_prob_norm = step_prob/torch.sum(step_prob)
# print(test_cat_exp)
# print(torch.sum(test_cat_exp))
# plt.imshow(test_cat_exp.detach().numpy()[:,:])
# plt.colorbar()
# plt.show()
# Flatten the probability surface
flat_step_prob_norm = step_prob_norm.flatten().detach().numpy()
# print(flat_prob_surface)
# Generate the corresponding indices for the flattened array
indices = np.arange(flat_step_prob_norm.size)
# Sample from the flattened probability surface
sampled_index = np.random.choice(indices, p=flat_step_prob_norm)
# Convert the sampled index back to 2D coordinates
sampled_coordinates = np.unravel_index(sampled_index, step_prob_norm.shape)
# THE Y COORDINATE COMES FIRST in the sampled coordinates
new_px = origin_xs[0] + sampled_coordinates[1]
new_py = origin_ys[0] + sampled_coordinates[0]
# Convert geographic coordinates to pixel coordinates
new_x, new_y = raster_transform * (new_px, new_py)
# Place the sampled location at a random point within the cell (rather than the centre)
# new_x = new_x + np.random.uniform(-12.5, 12.5)
# new_y = new_y + np.random.uniform(-12.5, 12.5)
# instead of the uniform, sample from a normal distribution with mean 0 and sd 6.5,
# which are the parameters where the cell contains ~ 95% of density
# if it's outside the bounds of the cell, resample
while True:
jitter_x = np.random.normal(0, 6.5)
if -12.5 <= jitter_x <= 12.5:
break
# Sample jitter for new_y and ensure it is within bounds
while True:
jitter_y = np.random.normal(0, 6.5)
if -12.5 <= jitter_y <= 12.5:
break
# Add the valid jitter to new_x and new_y
new_x = new_x + jitter_x
new_y = new_y + jitter_y
# print(new_x, new_y)
# return new_x, new_y, step_prob
return new_x, new_y, hab_log_prob, move_log_prob, step_log_prob, sampled_coordinates[1], sampled_coordinates[0]
Call the function¶
In [60]:
global_raster_list = [ndvi_global_norm[which_ndvi,:,:], canopy_global_norm, herby_global_norm, slope_global_norm]
x2 = x2_full[3,:].unsqueeze(0)
# x3 = x2
bearing_step = bearing[3].unsqueeze(0).unsqueeze(0)
print(x2)
print(bearing.shape)
simulate_trajectory(global_raster_tensors=global_raster_list,
scalars_to_grid=x2,
# additional_data_tensor=x3,
bearing=bearing_step,
window_size=101,
x_loc=41969,
y_loc=-1.435671e+06,
global_raster_transform=raster_transform)
tensor([[ 0.7071, 0.7071, -0.9942, 0.1074]]) torch.Size([1000])
Out[60]:
(41294.6246129592,
-1435876.496104774,
tensor([[-10.9798, -10.8766, -10.6232, ..., -10.8568, -11.0595, -11.0156],
[-11.2253, -10.9846, -10.5650, ..., -10.8998, -11.1774, -10.9955],
[-10.9037, -10.5621, -9.9899, ..., -10.4822, -10.9105, -10.8701],
...,
[-10.6507, -10.2017, -9.7033, ..., -9.8946, -10.2207, -10.4745],
[-11.0968, -10.9636, -10.6533, ..., -10.7428, -10.8820, -10.8480],
[-11.2361, -11.2483, -11.1158, ..., -11.1469, -11.1915, -11.0522]],
dtype=torch.float64, grad_fn=<SqueezeBackward0>),
tensor([[-15.7534, -15.6741, -15.5955, ..., -15.2336, -15.3086, -15.3843],
[-15.6785, -15.5983, -15.5188, ..., -15.1534, -15.2293, -15.3059],
[-15.6043, -15.5233, -15.4429, ..., -15.0738, -15.1506, -15.2282],
...,
[-15.6471, -15.5668, -15.4871, ..., -15.1587, -15.2349, -15.3118],
[-15.7220, -15.6425, -15.5637, ..., -15.2389, -15.3142, -15.3902],
[-15.7976, -15.7190, -15.6410, ..., -15.3198, -15.3942, -15.4692]],
dtype=torch.float64, grad_fn=<SqueezeBackward0>),
tensor([[-26.7332, -26.5507, -26.2187, ..., -26.0904, -26.3681, -26.3999],
[-26.9038, -26.5830, -26.0838, ..., -26.0531, -26.4066, -26.3014],
[-26.5080, -26.0854, -25.4327, ..., -25.5560, -26.0611, -26.0983],
...,
[-26.2979, -25.7685, -25.1904, ..., -25.0533, -25.4555, -25.7863],
[-26.8188, -26.6061, -26.2170, ..., -25.9818, -26.1962, -26.2381],
[-27.0337, -26.9673, -26.7567, ..., -26.4667, -26.5857, -26.5215]],
dtype=torch.float64, grad_fn=<SqueezeBackward0>),
23,
58)
Generate trajectory¶
Setup parameters¶
In [61]:
# Get today's date for saving figures
today_date = datetime.today().strftime('%Y-%m-%d')
# Setup the simulation parameters
n_steps = 3000
starting_yday = 206
# starting location of buffalo 2005
start_x = 41969.310875
start_y = -1.435671e+06
global_raster_list = [ndvi_global_norm, canopy_global_norm, herby_global_norm, slope_global_norm]
print(global_raster_list[0].shape)
window_size = 101
global_raster_transform = raster_transform
torch.Size([24, 2280, 2400])
Create simulation inputs from the parameters¶
In [62]:
x = np.repeat(0., n_steps)
y = np.repeat(0., n_steps)
# print(x)
x[0], y[0] = start_x, start_y
# print(x, y)
# Create sequence of steps
step = range(1, n_steps)
# print(step)
# hour of the day (hour) sequence
hour_t2 = np.resize(range(24), n_steps)
# print(hour_t2)
# convert to sine and cosine
hour_t2_sin = np.sin(2*np.pi*hour_t2/24)
hour_t2_cos = np.cos(2*np.pi*hour_t2/24)
# Create the day of the year sequences
yday_t2 = np.repeat(range(starting_yday, starting_yday + 365), 24)
yday_t2 = np.resize(yday_t2, n_steps)
# print(yday_t2)
# convert to sine and cosine
yday_t2_sin = np.sin(2*np.pi*yday_t2/365)
yday_t2_cos = np.cos(2*np.pi*yday_t2/365)
# bearing vector
bearing = np.repeat(0., n_steps)
# Convert lists to PyTorch tensors
hour_t2_tensor = torch.tensor(hour_t2).float()
hour_t2_sin_tensor = torch.tensor(hour_t2_sin).float()
hour_t2_cos_tensor = torch.tensor(hour_t2_cos).float()
yday_t2_tensor = torch.tensor(yday_t2).float()
yday_t2_sin_tensor = torch.tensor(yday_t2_sin).float()
yday_t2_cos_tensor = torch.tensor(yday_t2_cos).float()
bearing_tensor = torch.tensor(bearing).float()
# Stack tensors column-wise to create a tensor of shape [n_steps, 4]
x2_full = torch.stack((hour_t2_sin_tensor, hour_t2_cos_tensor, yday_t2_sin_tensor, yday_t2_cos_tensor), dim=1)
# x3_full = x2_full
# print(head(x2_full))
# List to hold filenames of saved images
filenames_hab = []
filenames_move = []
filenames_step = []
Select the appropriate NDVI layer from the yday¶
In [63]:
# Initialize variables to cache the previous yday and month index
previous_yday = None
month_index = None
day_to_month_index(500)
yday_sequence_1 = np.repeat(range(starting_yday, starting_yday + 365), 24)
print(yday_sequence_1)
print(len(yday_sequence_1))
[206 206 206 ... 570 570 570] 8760
Trajectory loop¶
In [64]:
for i in range(1, n_steps):
x_loc = x[i-1]
y_loc = y[i-1]
# calculate the bearing from the previous location
if i > 1:
bearing_rad = np.arctan2(y[i-1] - y[i-2], x[i-1] - x[i-2])
else:
bearing_rad = np.random.uniform(-np.pi, np.pi)
# Store the bearing in the vector
bearing[i-1] = bearing_rad
# print("Bearing[i-1]", bearing[i-1])
bearing_tensor = torch.tensor(bearing[i-1]).unsqueeze(0).unsqueeze(0)
x2 = x2_full[i-1,:].unsqueeze(dim=0)
# print(x2)
# Determine the month index based on the day of the year
day_of_year = yday_t2[i-1]
if day_of_year != previous_yday:
month_index = day_to_month_index(day_of_year)
previous_yday = day_of_year
# print(f'Day of the year: {day_of_year}')
month_index = day_to_month_index(day_of_year)
# print(f'Month index: {month_index}')
global_raster_list = [ndvi_global_norm[month_index,:,:], canopy_global_norm, herby_global_norm, slope_global_norm]
sim_outputs = simulate_trajectory(global_raster_tensors=global_raster_list,
scalars_to_grid=x2,
# additional_data_tensor=x2,
bearing=bearing_tensor,
window_size=101,
x_loc=x_loc,
y_loc=y_loc,
global_raster_transform=global_raster_transform)
new_x, new_y, hab_log_prob, move_log_prob, step_log_prob, px, py = sim_outputs
# print(new_x, new_y)
# print(px, py)
new_bearing = np.arctan2(new_y - y_loc, new_x - x_loc)
# print(new_bearing)
x[i] = new_x
y[i] = new_y
# ### PLOTTING
# # Create the mask for x and y coordinates
# x_mask = np.ones_like(hab_density)
# y_mask = np.ones_like(hab_density)
# # mask out cells on the edges that affect the colour scale
# x_mask[:, :3] = -np.inf
# x_mask[:, 98:] = -np.inf
# y_mask[:3, :] = -np.inf
# y_mask[98:, :] = -np.inf
# ## Habitat probability surface
# # normalise the raster for plotting
# hab_log_prob_norm = hab_log_prob/-torch.sum(hab_log_prob)
# # convert to numpy array
# hab_log_prob_norm = hab_log_prob_norm.detach().numpy()[:,:]
# # Create the mask for x and y coordinates
# hab_log_prob_norm = hab_log_prob_norm * x_mask * y_mask
# # Save the figure
# filename_hab = f"outputs/dl_prob_maps/ta_mix/id{buffalo_id}_hab_log_prob_norm_{today_date}_{i}.png"
# plt.figure() # Create a new figure
# plt.imshow(hab_log_prob_norm)
# plt.colorbar()
# plt.draw() # Ensure the plot is rendered
# plt.savefig(filename_hab, dpi=300)
# # plt.show()
# plt.close() # Close the figure to free memory
# # Add filename to the list
# filenames_hab.append(filename_hab)
# ### Movement probability surface
# # normalise the raster for plotting
# move_log_prob_norm = move_log_prob/-torch.sum(move_log_prob)
# # convert to numpy array
# move_log_prob_norm = move_log_prob_norm.detach().numpy()[:,:]
# # Save the figure
# filename_move = f"outputs/dl_prob_maps/id{buffalo_id}_move_log_prob_norm_{today_date}_{i}.png"
# plt.figure() # Create a new figure
# plt.imshow(move_log_prob_norm)
# plt.colorbar()
# plt.draw() # Ensure the plot is rendered
# plt.savefig(filename_move, dpi=300)
# # plt.show()
# plt.close() # Close the figure to free memory
# # Add filename to the list
# filenames_move.append(filename_move)
# ### Step selection probability surface
# # normalise the raster for plotting
# log_prob_surface_norm = step_log_prob/-torch.sum(step_log_prob)
# # convert to numpy array
# log_prob_surface_norm = log_prob_surface_norm.detach().numpy()[:,:]
# # Set the pixel of the next step, which is at (x, y) to 0
# log_prob_surface_norm[py, px] = -0.0001
# # Create the mask for x and y coordinates
# log_prob_surface_norm = log_prob_surface_norm * x_mask * y_mask
# # Save the figure
# filename = f"outputs/dl_prob_maps/ta_mix/id{buffalo_id}_step_log_prob_norm_{today_date}_{i}.png"
# plt.figure() # Create a new figure
# plt.imshow(log_prob_surface_norm)
# plt.colorbar()
# plt.draw() # Ensure the plot is rendered
# plt.savefig(filename, dpi=300)
# # plt.show()
# plt.close() # Close the figure to free memory
# # Add filename to the list
# filenames_step.append(filename)
# ### Step selection probability surface (exp)
# # normalise the raster for plotting
# prob_surface_norm = torch.exp(step_log_prob)/torch.sum(torch.exp(step_log_prob))
# # log_prob_surface_norm = step_log_prob/-torch.sum(step_log_prob)
# # convert to numpy array
# prob_surface_norm = prob_surface_norm.detach().numpy()[:,:]
# # Set the pixel of the next step, which is at (x, y) to 0
# prob_surface_norm[py, px] = 0.0
# # Save the figure
# filename = f"outputs/dl_prob_maps/ta_mix/id{buffalo_id}_prob_surface_{today_date}_{i}.png"
# plt.figure() # Create a new figure
# plt.imshow(prob_surface_norm)
# plt.colorbar()
# plt.draw() # Ensure the plot is rendered
# # plt.savefig(filename, dpi=300)
# plt.show()
# plt.close() # Close the figure to free memory
# # Add filename to the list
# filenames_step.append(filename)
Checking the model outputs¶
In [65]:
# print(bearing)
plt.hist(bearing[x>0], bins=50) #, edgecolor='black'
plt.title('Bearing values')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.show()
print(bearing.shape)
# Calculate the turning angle (difference between consecutive values)
ta = np.diff(bearing[x>0])
ta = np.insert(ta, 0, 0)
plt.hist(ta, bins=50) #, edgecolor='black'
plt.title('Turning angle values')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.show()
ta_corr = np.where(ta > np.pi, ta - (2 * np.pi),
np.where(ta < -np.pi, ta + (2 * np.pi), ta))
plt.hist(ta_corr, bins=50) #, edgecolor='black'
plt.title('Turning angle values')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.show()
print(ta.shape)
step_log_prob_norm = step_log_prob/torch.sum(step_log_prob)
# print(step_log_prob_norm)
plt.imshow(step_log_prob.detach().numpy()[:,:])
plt.colorbar()
plt.show()
# check for NaN values
# print(torch.isnan(step_log_prob_norm).any())
step_prob = torch.exp(step_log_prob)
step_prob_norm = step_prob/torch.sum(step_prob)
# print(torch.sum(step_prob_norm))
plt.imshow(step_prob_norm.detach().numpy()[:,:])
plt.colorbar()
plt.show()
# # Create a histogram
# plt.hist(step_prob_norm.detach().numpy()[:,:], bins=100) #, edgecolor='black'
# # Add title and labels
# plt.title('Histogram Example')
# plt.xlabel('Value')
# plt.ylabel('Frequency')
# Show the plot
plt.show()
# ndvi_global_norm = (ndvi_global_tens - ndvi_min) / (ndvi_max - ndvi_min)
(3000,)
(3000,)
Plot the simulated trajectory¶
In [66]:
# Actual number of locations before there were NaNs
# print(x)
print(x[x>0].shape[0])
# Create a figure and axis with matplotlib
fig, ax = plt.subplots(figsize=(7.5, 7.5))
# Plot the raster
show(ndvi_global[which_ndvi,:,:], transform=raster_transform, ax=ax, cmap='viridis')
# Set the title and labels
ax.set_title('NDVI')
ax.set_xlabel('Easting')
ax.set_ylabel('Northing')
# Plot the simulated trajectory
no_sim_points = np.min([x[x>0].shape[0], y[y<0].shape[0]])
print(no_sim_points)
plt.plot(x[1:no_sim_points], y[1:no_sim_points], color = 'red')
plt.plot(x[1:no_sim_points], y[1:no_sim_points], color = 'red')
# Show the plot
plt.show()
3000 3000
Write the trajectory to a csv¶
Only run this once otherwise it will create duplicates¶
In [67]:
# Combine vectors into a DataFrame
trajectory_df = pd.DataFrame({'x': x[x>0],
'y': y[x>0],
'hour': hour_t2[x>0],
'yday': yday_t2[x>0],
'bearing': bearing[x>0],
'ta': ta})
n_steps_actual = x[x>0].shape[0]
# Save the DataFrame to a CSV file
index = 1
csv_filename = f'outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id{buffalo_id}_{n_steps_actual}steps_{index}_{today_date}.csv'
# Check if the file already exists and find a new name if necessary
while os.path.exists(csv_filename):
csv_filename = f'outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id{buffalo_id}_{n_steps_actual}steps_{index}_{today_date}.csv'
index += 1
print(csv_filename)
trajectory_df.to_csv(csv_filename, index=True)
outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_3000steps_5_2024-09-29.csv
Multiple trajectories in a loop¶
In [68]:
# Setup the simulation parameters
# here we just run a couple simulations to show that it's working
n_trajectories = 5
n_steps = 1000
starting_yday = 206
# global_raster_list = [ndvi_global_norm, canopy_global_norm, herby_global_norm, slope_global_norm]
window_size = 101
global_raster_transform = raster_transform
# Looping over individuals
for j in range(1, n_trajectories+1):
# Setup the vectors used in the prediction function
x = np.repeat(0., n_steps)
y = np.repeat(0., n_steps)
x[0], y[0] = start_x, start_y
# Create sequence of steps
step = range(1, n_steps)
# hour of the day (hour) sequence
hour_t2 = np.resize(range(24), n_steps)
# convert to sine and cosine
hour_t2_sin = np.sin(2*np.pi*hour_t2/24)
hour_t2_cos = np.cos(2*np.pi*hour_t2/24)
# Create the day of the year sequences
yday_t2 = np.repeat(range(starting_yday, starting_yday + 365), 24)
yday_t2 = np.resize(yday_t2, n_steps)
# convert to sine and cosine
yday_t2_sin = np.sin(2*np.pi*yday_t2/365)
yday_t2_cos = np.cos(2*np.pi*yday_t2/365)
# bearing vector
bearing = np.repeat(0., n_steps)
# Convert lists to PyTorch tensors
hour_t2_tensor = torch.tensor(hour_t2).float()
hour_t2_sin_tensor = torch.tensor(hour_t2_sin).float()
hour_t2_cos_tensor = torch.tensor(hour_t2_cos).float()
yday_t2_tensor = torch.tensor(yday_t2).float()
yday_t2_sin_tensor = torch.tensor(yday_t2_sin).float()
yday_t2_cos_tensor = torch.tensor(yday_t2_cos).float()
bearing_tensor = torch.tensor(bearing).float()
# Stack tensors column-wise to create a tensor of shape [n_steps, 4]
x2_full = torch.stack((hour_t2_sin_tensor, hour_t2_cos_tensor, yday_t2_sin_tensor, yday_t2_cos_tensor), dim=1)
# Initialize variables to cache the previous yday and month index
previous_yday = None
month_index = None
# simulation loop
for i in range(1, n_steps):
x_loc = x[i-1]
y_loc = y[i-1]
# calculate the bearing from the previous location
if i > 1:
bearing_rad = np.arctan2(y[i-1] - y[i-2], x[i-1] - x[i-2])
else:
bearing_rad = np.random.uniform(-np.pi, np.pi)
# store bearing in the vector
bearing[i-1] = bearing_rad
bearing_tensor = torch.tensor(bearing[i-1]).unsqueeze(0).unsqueeze(0)
x2 = x2_full[i,:].unsqueeze(dim=0)
# Determine the month index based on the day of the year
day_of_year = yday_t2[i-1]
if day_of_year != previous_yday:
month_index = day_to_month_index(day_of_year)
previous_yday = day_of_year
# print(f'Day of the year: {day_of_year}')
month_index = day_to_month_index(day_of_year)
# print(f'Month index: {month_index}')
global_raster_list = [ndvi_global_norm[month_index,:,:], canopy_global_norm, herby_global_norm, slope_global_norm]
sim_outputs = simulate_trajectory(global_raster_tensors=global_raster_list,
scalars_to_grid=x2,
bearing=bearing_tensor,
window_size=101,
x_loc=x_loc,
y_loc=y_loc,
global_raster_transform=global_raster_transform)
new_x, new_y, hab_log_prob, move_log_prob, step_log_prob, px, py = sim_outputs
x[i] = new_x
y[i] = new_y
# save the data frames individually
# Combine vectors into a DataFrame
trajectory_df = pd.DataFrame({'x': x[x>0],
'y': y[x>0],
'hour': hour_t2[x>0],
'yday': yday_t2[x>0],
'bearing': bearing[x>0]})
n_steps_actual = x[x>0].shape[0]
# Save the DataFrame to a CSV file
index = j
csv_filename = f'outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id{buffalo_id}_{n_steps_actual}steps_{index}_{today_date}.csv'
# Check if the file already exists and find a new name if necessary
while os.path.exists(csv_filename):
csv_filename = f'outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id{buffalo_id}_{n_steps_actual}steps_{index}_{today_date}.csv'
index += 1
print(csv_filename)
trajectory_df.to_csv(csv_filename, index=True)
outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_1000steps_6_2024-09-29.csv outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_1000steps_7_2024-09-29.csv outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_1000steps_8_2024-09-29.csv outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_1000steps_9_2024-09-29.csv outputs/dl_trajectories/ta_mix_monthly_ndvi/dl_trajectory_TAmix_id2005_1000steps_10_2024-09-29.csv